
function res = trebufun(theta,thetadot,phi,phidot,t,L1, L2,M1,M2,R,g)
phi = (phi/180)*pi;
theta = (theta/180) *pi;
trebvars = [theta;phi;thetadot;phidot] ;
[T,K] = ode45(@(t,trebvars)acceleration(trebvars,L1,L2,M1,M2,R,g),[0,t],trebvars);
res = [T,K];
end

function acc=acceleration(trebvars,L1,L2,M1,M2,R,g)
 theta = trebvars(1);
phi = trebvars(2);
thetadot = trebvars(3);
phidot = trebvars(4);
leftover = [L1*L2*M2*(thetadot^2)*sin(theta-phi)+g*M2*L2*sin(phi)
            -L1*L2*M2*(phidot^2)*sin(theta-phi)+g*R*M1*cos(theta)+g*L1*M2*sin(theta)]; %right side of matrix equation, top is phi
        
dotstuff = [L1*L2*M2*cos(theta-phi),-M2*L2
           -(M1*R^2+M2*L1^2),L1*L2*M2*cos(theta-phi)] ;                %left side of matrix equation, top is phi
      
doubledot =  dotstuff\leftover  ;              %calulates phi and theta umblat, in that order
    

acc = [thetadot;phidot;doubledot];
end